#ifndef U_H
#define U_H

#include <cblas.h> // 对于BLAS
#include <lapacke.h> // 对于LAPACK
#include <vector>
#include <cmath>
#include <iostream>
#include <string>

std::vector<double> copyToa(double b[], int n) {
    return std::vector<double>(b, b + n);
}

double B_u(int i, int degree, double t, std::vector<double>& points) {
    if (degree == 0) {
        if (points[i] <= t && t < points[i + 1]) return 1.0;
        return 0.0;
    }
    double left = (t - points[i]) / (points[i + degree] - points[i]);
    double right = (points[i + degree + 1] - t) / (points[i + degree + 1] - points[i + 1]);
    return left * B_u(i, degree - 1, t, points) + right * B_u(i + 1, degree - 1, t, points);
}

double USplineReturn(double t,  std::vector<double>& x,  std::vector<double>& y, double fa, double fb) {
    int degree = 2;
    int N = x.size() + 1;
    double L = x[N - 2] - x[0];

    // 扩展 y 向量
    std::vector<double> y1 = y;
    y1.insert(y1.begin(), fa);
    y1.push_back(fb);

    // 定义节点向量 p
    std::vector<double> p = x;
    p.insert(p.begin(), x[0] - L / (N - 2) / 2);
    p.push_back(x[N - 2] + L / (N - 2) / 2);

    // 定义插值点
    std::vector<double> points = x;
    for (int i = 0; i < degree; ++i) {
        points.insert(points.begin(), x[0] - (i + 1) * L / (N - 2));
        points.push_back(x[N - 2] + (i + 1) * L / (N - 2));
    }

    int n = N + 1;  // 需要解的方程数目

    // 确保矩阵 A 和 b 的大小正确
    double A[n * n];  // 存储系数矩阵 A，大小为 n*n
    double b[n];   // 存储常数向量 b，大小为 n
    std::fill(A, A + n * n, 0.0);
    std::fill(b, b + n, 0.0);

    int ipiv[n];  // 主元素的排列
    int lda = n;  // A 的 leading dimension
    int ldb = 1;  // b 的 leading dimension
    int info;     // LAPACKE 返回的状态

    // 填充矩阵 A
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            A[i * n + j] = B_u(j, degree, p[i], points);
        }
    }

    // 填充常数向量 b
    for (int i = 0; i < n; ++i) {
        b[i] = y1[i];
    }

    /*// 打印矩阵 A
    std::cout << "Matrix A:" << std::endl;
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            std::cout << A[i * n + j] << " ";
        }
        std::cout << std::endl;
    }

    // 打印常数向量 b
    std::cout << "Solution to the system (first derivatives):" << std::endl;
    for (int i = 0; i < n; ++i) {
        std::cout << b[i] << std::endl;
    }*/

    // 使用 LAPACKE 求解线性方程组 Ax = b
    info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, 1, A, lda, ipiv, b, ldb);

    // 检查求解结果
    if (info == 0) {
        std::vector<double> a = copyToa(b, n);

        double result = 0.0;
        // 计算加权基函数的结果
        for (int i = 0; i < n; ++i) {
            result += a[i] * B_u(i, degree, t, points);
        }
        return result;
    } else {
        std::cerr << "Error in solving the system, info = " << info << std::endl;
    }

    return 0.0;
}

#endif
